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Abstract 

In this work we develop a complete variational many-body theory for a system of N trapped 
bosons interacting via a general two-body potential. The many-body solution of this system is 
expanded over orthogonal many-body basis functions (configurations). In this theory both the 
many-body basis functions and the respective expansion coefiicients are treated as variational 
parameters. The optimal variational parameters are obtained self-consistently by solving a coupled 
system of non-eigenvalue - generally integro-differential - equations to get the one-particle functions 
and by diagonalizing the secular matrix problem to find the expansion coefficients. We call this 
theory multi-configurational Hartree for bosons or MCIIB(M), where M specifies explicitly the 
number of one-particle functions used to construct the configurations. General rules for evaluating 
the matrix elements of one- and two-particle operators are derived and applied to construct the 
secular Hamiltonian matrix. We discuss properties of the derived equations. We show that in 
the limiting cases of one configuration the theory boils down to the well-known Gross-Pitaevskii 
and the recently developed multi-orbital mean-fields. The invariance of the complete solution with 
respect to unitary transformations of the one-particle functions is utilized to find the solution with 
the minimal number of contributing configurations. 

In the second part of our work we implement and apply the developed theory. It is demonstrated 
that for any practical computation where the configurational space is restricted, the description of 
trapped bosonic systems strongly depends on the choice of the many-body basis set used, i.e., self- 
consistency is of great relevance. As illustrative examples we consider bosonic systems trapped in 
one- and two-dimensional symmetric and asymmetric double-well potentials. We demonstrate that 
self-consistency has great impact on the predicted physical properties of the ground and excited 
states and show that the lack of self-consistency may lead to physically wrong predictions. The 
convergence of the general MCHB(M) scheme with a growing number M is validated in a specific 
case of two bosons trapped in a symmetric double-well. 

PACS numbers: 03.75.Hh,03.65.Ge,03.75.Nt 
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I. INTRODUCTION 



The first experimental realizations of Bose-Einstein condensation in trapped ultra-cold 
atomic clouds renewed a great interest in the experimental and theoretical de- 

scriptions of this phenomenon. Modern experimental setups utilize magnetic 0|, electric 

and optical dipole fields or their combinations to control the trapping and guiding 
of the ultra-cold atoms. The number of condensed atoms in these experiments varies from 
several dozens l7f to several millions Magnetically-tunable Feshbach resonances 



make it possible to control the strength and sign of the inter-particle interactions. All 
these experimental tools may be used to design a bosonic system jo], |lO|, and to study its 
time-independent and time-dependent properties. 

From the theoretical point of view we have to study the properties of the collection of 
interacting many-electron atoms immersed in a time-dependent crossing electric and mag- 
netic fields. This very complicated quan tum mechanical many-body problem is replaced 
usually by a model Hamiltonian jl^. Il .4 ll^ Typically, the diluteness of the atomic 
cloud allows to consider the atoms as point-like particles with pairwise interaction and to 
neglect three-body and higher order collisions. The crossing electric and magnetic fields 
are replaced by an effective external trap potential. Within these assumptions the original 
system is modelled as a collection of massive point-like particles interacting via repulsive 
or attractive pairwise inter-particle interaction immersed in an external trap potential of 
known geometry. 

However, there are only a few examples of the model many-body Hamiltonians the exact 
solutions of which are known, see, e.g., Refs. 0,0,0 and references therein. For general 
inter-particle interactions and trap geometries the problem can be attacked only within the 
framework of approximate and numerical methods. The most popular one is the variational 
approach, where the form or a specific ansatz of the trial many-body function is postulated. 
This ansatz depends on several parameters, and to solve the problem means to find the op- 
timal set of the parameters which would minimize the expectation value of the Hamiltonian. 
Clearly, the more parameters are involved the closer the obtained solution to the exact true 
many-body function is. 

One of the widely and successfully used approximations is the Gross-Pitaevskii (GP) 
ansatz Q], where it is assumed that each boson resides in a single spatial function, i.e.. 
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the many-body bosonic solution is presented as a product of identical one-particle functions 
(GP orbital): 

^'(fi,f2, . . . ,fAr) = ip{fi)ip{f2) ■ ■■ip{rN). (1) 

This ansatz has only one variational parameter - the shape of the GP orbital y?. Orbitals 
with different shapes give different approximations to the true many-body solution. The 
"optimal" orbital is obtained by minimizing the expectation value of the Hamiltonian with 
respect to which is equivalent to solving the very known GP equation. The GP solution 
is self- consistent, i.e., the shape of the GP orbital depends on the number of bosons and 
strength of the inter-particle interaction in addition to the geometry of the external t rap 
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potential. However, despite the great success of the GP ansatz, see Refs.j 
and references therein, there are ample physical situations this one-orbital mean-field theory 
cannot describe, such as the depletion and fragmentation of trapped condensates and ap- 
pearance of Mott-insulator phases of cold bosonic atoms in optical lattices. The natural way 
to resolve these difficulties is to take a more general many-body ansatz and to go, therefore, 
beyond Gross-Pitaevskii theory. 
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Recently, a more general mean-field theory for bosons has been put forward 
The multi-orbital, or best mean-field (BMP) theory as it is also called has been derived by 
considering the many-body ansatz where ni bosons reside in one orbital 0i, n2 bosons in a 
second orbital 02 and so on: 

■^{ri, ...,rN)= 501 (ri) ■ ■ ■ 0l(r;j02(^ni+l) ■ ■ ■ 02(^^1+^2) " ■ ■ (/"Af (^ni+n2+-+n„), (2) 

where S is the symmetrizing operator required for bosons. In the multi-orbital mean-field 
theory the shapes of the one-particle functions 0j, the occupation numbers as well as 
the number M of one-particle functions are considered as variational parameters. The opti- 
mal parameters of this many-body function are obtained by solving a respective system of 
coupled non-linear equations. The resulting mean-field solution is also self- consistent, and 
depends on the number of bosons and strength of the inter-particle interaction in addition 
to the geometry of the external trap potential. The BMP theory was used, among several 
applications and predictions, to show 2^ that with increasing inter-particle interaction, the 
ground state of a trapped bosonic system on its pathway from condensation towards fermion- 
ization, gradually passes through two-, three- and up to M-fold stages of the fragmentation. 

n 

The multi-orbital mean-field theory, when applied to the optical lattices, predicts |2^] the 



transition from the super-fluid to the Mott-insulator phase and reveals a variety of new 
Mott-Insulator phases. 

The self-consistent GP and BMP theories considered above successfully describe the main 
features of the condensation, fragmentation and fermionization phenomena. However, the 
mean-field ansatzes used are made of only one many-body function of known type (a per- 
manent) and are, therefore, incapable - by construction - to describe the depletion and 
fluctuations of the many-body states. To improve the many-body description further it is 
natural to go beyond mean-fields and take as an ansatz a linear combination of several known 
N-body basis functions: 



where is a set of N-body basis functions and {Ci} are the expansion coefficients. In 

quantum chemistry, Eq.Q is called configuration interaction (CI) expansion [2^, because 
every many-body basis function (configuration) is attributed to a simple physical situation. 
Por bosonic systems, a configuration may represent a condensate, i.e., the state where all 
bosons reside in the same orbital, an excited state where N — i bosons remain condensed and 
i bosons are excited out of a condensate, a two-fold fragmented state where a macroscopic 
number of bosons ni reside in one orbital and n2 = N — ni bosons in another, and so on. The 
variational problem of finding the unknown expansion coefficients in this case is reduced 
to the diagonalization of the respective secular Hamiltonian matrix [2^. This method is 
also known as exact diagonalization technique, because if all possible configurations are 
considered, i.e., the N-body basis is complete, then this expansion is exact, irrespective 
to the particular choice of the many-body basis functions However, the secular 

Hamiltonian matrix in this case is an infinite matrix. 

To make real computations tractable, the number of configurations in the expansion 
Eq.Q, of course, has to be truncated. The configurational space in this case spans a re- 
stricted subspace of the Hilbert space. The exact diagonalization studies within truncated 
configurational spaces have been used to pro vide a many-body description of repulsive and 
attractive condensates, see, e.g., 

Refs.QBQ- 

In these investigations, many-body basis 
functions comprised of a priori fixed orbitals have been utilized to study the properties of the 
bosonic system as a function of the inter-particle interaction strength. However, it is clear 
that exact diagonalizations performed within different truncations of the CI expansion will 




(3) 
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give different results. Moreover, truncated CI expansions of the same size but utilizing dif- 
ferent sets of orbitals will also give different results in general. Indeed, it was demonstrated 
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29| that the total energy obtained by the exact diagonalization in a restricted configura- 



tional space can sometimes be even worse than the self-consistent GP mean-field energy. In 
these cases we encounter the situation that due to the implemented self-consistency, a single 
mean-field function provides a better description than a very large fixed-orbital many-body 
CI expansion. 

A question addressed in this work is how the choice of the many-body basis functions 
impacts the results obtained within a restricted configurational space. By comparing the 
many-body results obtained within different basis sets we can find the energetically most 
favorable one, i.e., the best many-body basis set. The main goal of the present investigation 
is to formulate a many-body theory which would provide the best set of many-body basis 
functions in a desired, i.e., truncated, subspace of the Hilbert space. To achieve this goal, 
we apply a general variational principle where we treat both the set of many-body basis 
functions and the set of the expansion coefficients {Cj} appearing in the expansion 

Eq.(jni) as variational parameters which are to be optimized. This results, as we shall see 
below, in a many-body theory for bosonic systems with complete self-consistency, which 
we refer to as Multi-Configurational Hartree for Bosons or MCHB(M), where M specifies 
explicitly the number of one-particle functions used to construct the configurations. 

Some hints that self-consistency is useful and important in attacking the many-boson 
problem beyond mean-field have already been addressed in the literature in the context 



of symmetric doub 
Spekkens and Sipe 



e-well potentials and two-orbital configuration-interaction expansions. 



30[ provide an approximative analytic as well as numerical solutions for 



the bosonic system trapped in symmetric double-well potentials in the regime where the 
inter-particle interaction can be treated perturbatively. More recently, Reinhard 31| and 
coworkers have combined a partial nonlinear optimization of the many-body basis functions 
with a linear variational principle for the expansion coefficients to describe ground and 
excited states of bosons trapped in a symmetric double-well. Here, we would like to stress 
that the many-body variational theory which we develop in this paper is complete, because 
we fully optimize all expansion coefficients as well as all the many-body basis functions 
appearing in the expansion Eq.©, and (/enera/ because it is valid for any trap geometry, for 
any physical shape of the inter-particle interaction and in any dimension. 
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The structure of the paper is as follows. In Sec III Al a trapped N-bosonic system inter- 
acting via a general pairwise potential is considered. The multi-configurational expansion 
Eq.Q is used as an ansatz for the true many-body wave function of this system. The 
variational principle is used in Sec lllBl to formulate the general equations which allow one 
to find the best many-body basis functions and the corresponding expansion coefficients 
self-consistently in the desired active space. Sec JIIB H is devoted to the problem of finding 
the expansion coefficients. In this section we provide the general rules for evaluating matrix 
elements of one- and two-body operators between two general basis functions (permanents) 
and apply them to construct the secular Hamiltonian matrix. The problem of finding the op- 
timal basis functions is considered in Sec JllB2l where the working equations for the general 
case are derived in closed form using the elements of the reduced one- and two-body density 
matrices. Properties of the resulting equations are discussed in Sec JII CI In particular, we 
demonstrate that in the limiting cases of a single permanent, the derived equations boil down 
to the one-orbital Gross-Pitaevskii and multi-orbital mean-fields equations. Sec lIIIAl opens 
the second part of our work and provides explicitly the working equations of MCHB(2), i.e., 
the case where the basis functions are obtained as all possible permutations of N bosons over 
two orbitals. As an inter-particle potential we implement the popular contact interaction. 
The formalism is used to investigate the impact of self-consistency on many-body predictions. 
The ground state of N=1000 bosons trapped in a symmetric one-dimensional double-well 
trap is investigated in Sec lIII Bl and for an asymmetric one in Sec illl (^1 In Sec illlDl we 
show for these examples that the number of many-body basis functions contributing to the 
expansion in Eq.Q can be significantly reduced by appropriately choosing the orbitals used 
to construct the many-body basis functions. Excited states of the bosonic system trapped 
in a symmetric double-well trap are investigated in Sec lIII El The two-dimensional bosonic 
system trapped in symmetric double-well is studied in Sec lIII Fl The convergence of the 
MCHB(M) theory is addressed in Sec lIVI where we apply different levels M= 2, ■ ■ ■ , 10 of 
MCHB(M) theory to study ground state properties of two bosons trapped in a symmetric 
double-well trap. Finally, SecfVl summarizes our results and conclusions. 
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II. THEORY 



A. Preliminaries 

Consider a system of identical spinless bosons of mass m immersed in an external time- 
independent trap potential V{r) and interacting via a general pairwise interaction potential 
W{fi — fj) where is the position of the i-th particle. By using bosonic annihilation and 
creation operators bi and foj, — ^ij which are associated with a given set of orbitals 

{4>i}, the Hamiltonian of the system takes on the standard form in second quantization 
language: 

H = h + W 

h=Y^h,,b\b, , W=\Y. W^3kAb]hbk, (4) 

i,j = l i,j,k,l=l 

where the one- and two-body matrix elements read 

h^, = I <l^^i^)i-^^' + V{r))<P,{f)df (5) 

and 

Wi,M = j j <P:{^<l)*{f')W{f-f')Mf)M^)dfdf'. (6) 

Our general intention is to find time-independent solutions of this Hamiltonian in a 
form of expansion Q over basis functions. Every basis function being a many-body 
function must depend on the coordinates of all N bosons. The simplest way to con- 
struct it is to take a product of different orthogonal orbitals, so called Hartree product 
<l>j = r2, . . . , r]\f) = iS0i(ri)02(^2) ■ ■ ■ (pNi^jy), and to apply the symmetrizing operator 

S to fulfil the Bose statistic. When all the orbitals are identical, we obtain a GP-like many- 
body basis function (see Eq.(^) used to describe condensation. If some fraction of bosons 
resides in one orbital and the rest in another one, we deal with a basis function describing 
two-fold fragmentation, and so on. The general many-body basis function can be considered 
as one of the configurations resulting due to permutation of N bosons over M orbitals. In 
second quantization language this general many-body function, also known as permanent, 
reads 

Mri,r2, . . . , r-W) = , , ,\ =^{b\r (4)"' " " " (4)"*^ \vac) = |ni, ns, ng, ■ ■ ■ , Um) 

(7) 
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where M is a number of the one-particle functions, ni + n2 + + ■ ■ ■ + um = N ^ and \vac) 
is the vacuum. 

We recall that if the many-body basis set in expansion Q is complete, i.e., it spans the 

whole N-boson Hilbert space, then expansion Q is exact irrespective to the particular choice 

of the basis functions used. In the present formulation, this is achieved when the number 

of the one-particle functions M ^ oo. However, to make real computations tractable the 

size of the expansion (jS)) and hence the number of the one-particle basis functions M has to 

be restricted. Let us assume that the many-body basis set is formed by all possible 

configurations appearing as permutation of N bosons over M orbitals. The total number of 

/ M + N -l\ 

configurations, and therefore the size of the expansion in this case is . For 

example, for a system of N=1000 bosons and M=3 orbitals the total number of configurations 
is 501,501 while for M=4, the number of configurations is already 167,668,501, which would 
make practical computations impossible. Another consequence of the truncation is that 
different sets of orbitals used to construct many-body basis sets of the same size would lead 
to different approximations to the exact many-body \E'. As truncated CI expansions become 
very demanding with increasing N and/or M, it is of great advantage to exploit this property. 
Namely, it is desirable to find not only expansion coefficients but also the "best" set of one- 
particle basis functions. Combined together, these lead to the self-consistent optimal CI 
expansion. To fulfil this goal, we apply in this work the variational principle, which defines 
the energetically most favorable solution as the best one. 



B. The general variational approach 

Let us consider a system of N bosons and restrict the number of one-particle functions 
to M. Then, the trial many-body function (expansion Q) takes on the following form: 

1^) = X] C'(ni,n2,n3,-,nM)l^l''^2,%,--- = J^C'^ln), (8) 

ni,n2,--- ,njv/ n 

where ni, ?t,2, ■ ■ ■ , um are the number of bosons residing in orbitals 0i, 02, ■ ■ ■ 0Af- The sum- 
mation runs over all possible configurations n = (ni, n2, n^, ■ ■ ■ , hm), preserving total num- 
ber of particles ni+n2 + ■ ■ ■ + nM = N. The expectation value of the Hamiltonian evaluated 
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with this trial function reads: 



{mm = Y.ClC^,{n\H\n'). (9) 

n,n' 

This expression depends on two types of variational parameters: on the expansion coef- 
ficients {Cfi = C(ni,n2,n3,---,nM)} and on the particular choice of the one-particle functions 
01,02, ■ ■ -(pM = {(pi}- Our main goal is to find the values of these parameters for which 
(\l'|i^|\l/) is a minimum. The natural requirements on normalization of the trial many-body 
function (^E'l^E') = 1 and orthonormalization of all the one-particle functions {(j)i\(f)j) = 6ij 
allow us to formulate the minimization problem within Lagrange's method of undetermined 
multipliers. Here we have to stress that these are the only constrains applied. 
The Lagrange energy functional takes on the form: 

M M 

n 2=1 j=l 

where S appears due to normalization constrain of the trial many-body function (jHl), and 
fiij - due to orthonormalization of all one-particle functions. To find the minimum of 
this functional we put to zero all first derivatives of this functional with respect to every 
<^(rii,n2,n3,-,nM) appearing in the expansion Eq.®: 

dC[{CH},m] 



= J2Cn'{n\H\n') -SCri, Vn (11a) 



SC[{C^}m] ^ Q ^ _ ^^^10^) _ ^^,10,) . . . _ ^ = l,..., M.(llb) 

Here we used partial derivatives and functional derivatives to separate the variations with 
respect to the expansion coefficients and the one-particle functions. Such a separation is 
permitted because /ijj and (pi do not depend explicitly on C^. Eq. ljllap defines the expansion 
coefficients when the set of the one-particle functions is given, while Eq. ()llb|) finds the 
best, i.e., energetically most favorable one-particle functions when the set of the expansion 
coefficients is known. 

We can use the following strategy to solve the two-fold variational problem Eqs. ljllaj) and 
flllbj) self- consistently. Starting from some guess for the orbitals {(pi\ we construct the initial 
many-body basis set. Then we use these initial fixed-orbital many-body basis functions to 
build up the secular Hamiltonian matrix. As we shall see in the subsequent Sec JII B l| the 
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first part of tlie variational principle Eq. ()llap . i.e., the problem of finding the unknown 
coefficients C{ni,n2,n3,-- ,nM) "^^^ be reduced to the diagonalization of the secular Hamiltonian 
matrix. The expansion coefficients obtained as the components of the respective eigenvector 
are utilized in the second part of the variational procedure Eq. ()llb|) which provides as its 
output a new approximation to the one-particle functions {4>i}, see for details Sec ill B 2l This 
iterative scheme is repeated until convergence is achieved. We call the obtained optimal set 
of the one-particle functions and respective expansion coefficients the self-consistent solution. 

Since in this approach the many-body bosonic wave-function is presented as a sum of 
all possible configurations (symmetrized Hartree products) resulting as permutations of N 
bosons over M orbitals, we name this general variational method Multi-Configurational 
Hartree for Bosons or MCHB(M). Here M specifies explicitly the number of one-particle 
functions involved. 

1. Variations with respect to the expansion coefficients {C^} and general rules for evaluating 
matrix elements with permanents 

The quantities {n\H\n') appearing in Eq. ljllap are the matrix elements of a matrix Hamil- 
tonian Ti. This system of equations can be written in matrix notations as 

nC = SC, (12) 

where C is a column vector of the expansion coefficients C^. As usual, the problem of finding 
the expansion coefficients {Ca} for a given set of one-particle functions {(pi} boils down to 
an eigenvalue problem. This opens the possibility to attack not only the ground, but also 
excited states. 

The main issue now is to evaluate the matrix elements (n|if|n') with permanents. As 
introduced in Eq.(@]), the Hamiltonian is made of an one-particle operator h and of a two- 
body inter-particle interaction operator W. It is useful to treat these one-body {n\h\n') and 
two-body (fi|W^|n') interaction terms individually. 

In the following we report the general rules to evaluate the matrix elements of one-body 
and two-body operators. For simplicity these operators are denoted h and W although they 
need not be the constituents of the Hamiltonian. Let us distinguish between six generic 
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types of permanents: 



\Po) = 


1^1, 


n2, • • ■ 


,Tli,--- ,Tlj, 






1 J ''^Ji 










\P1) = 




^2, • • • 


,ni + l,--- 




-I,-- 


,nM) 




7ii + l;7ij - 


1;) 




1^2) = 




712, •• • 


, 7lj + 2, • • • 


rij 


_2,... 


,riM) 




7ij + 2; 7ij — 


2;) 




1^3) = 




n2, • • • 


, 'T'i + 2, • • • 


Uj 


-I,--- 


,nk - 


-1,- 


• • ,nM) = ; 


+ 2;nj - l;nfc - 


1;) 


1^4) = 




712, •• • 


,ni + l,--- 


Uj 


+ !,••• 


,nk - 


-2,- 


• • ,1^-1^) = ; 


7ij + l;nj + l;7ifc - 


2;) 


m = 




7l2, • • • 


,ni + l,--- 


Uj 


+ !,••• 


,rik - 


-1,- 


■■ ,ni-l,- 


• ,?T'm) 






|;n, 


+ l;rz^ 


j + l;nk- 1; 


ni 


-I;)- 










(13) 



The permanent |Pi) can be obtained from |Po) by excitation of a single boson from (f)j to 
(f)i. The other four permanents, IP2) to jPs), describe two-boson excitations of |Po). To 
shorten the notations we show only the occupation numbers of the involved orbitals. For 
instance, the configuration (; rii — 1; rij + 1; ) differs from n = (; nf, nj; ) by an excitation of 
a single boson from 0j to (f)j. We stress that only if two permanents differ by the excitation 
of at most two bosons, the Hamiltonian matrix element evaluated with these permanents is 
non-zero. All other permanents give vanishing matrix elements. 

The matrix elements of an one-body operator h are non-zero if the two permanents are 
either equal (diagonal contributions) or differ by an excitation of a single boson: 

M M 

{Po\h\Po) = {■,ni;nj;\ ^ Kpl>^J}i3\;ni]nj]) ^^haaUa, 

a,/3=l a=l 

(Pi|A|Po) = (;7ii + l;nj - 1; \h\;ni;nj;) = hij^/rii + ly^. (14) 
The diagonal matrix elements of a two-body operator W are: 

{Po\W\Po) = {■,ni;nj; \ ^ hVai3'ysbiblb^bs\;ni;nj;) 

a,l3,y,S 

^ M I ^ 

= ^^riairia- l)Waaaa + ^ nano{Waljafj + WajBIja) ■ (15) 

The matrix elements of a two-body operator W evaluated with permanents which differ by 
the excitation of one boson read: 

{Pi\W\Pq) = {■ni + l]nj-l-\W\]ni;nj-) 

I M 

= ^Vrk^^/n][ ^ naiWiaaj + W^iaja) + ^iM^mj + (% " ^Wijjj]. (16) 
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The permanents which differ by the excitations of two bosons give the following matrix 
elements of the two-body operator W: 

{P2\W\Po) = {■,n^ + 2;nj-2;\W\;ni;n,;) 

= \\l K- - l)^jV(«i + 2)(ni + l)Wiijj, 

{Pz\W\Pq) = {■,m + 2;nj - l;nk- l;\W\;ni;nj;) 

= a/ {rii + 2){ni + l)^ynjnkWiijk, 

{P4\W\Po) = {■,n, + l;nj + l;nk-2;\W\;nf,nf,) 

= V^k{nk - l)\J {rii + Vjijij + l)Wijkk, 

{P^\W\Po) = {■,ni + l;nj + l;nk - l;ni - l;\W\;ni;nj;) 

= ^ {m + l)(nj + l)^/^i{Wijki + Wjiki). (17) 

In summary, we have demonstrated that for any set of orthogonal one-particle functions 
{(pi} used to construct the many-body basis functions, i.e., the permanents, the variational 
problem of finding the unknown coefficients C{ni,n2,n3,---,nAi) is reduced to the diagonalization 
of the Hamiltonian secular matrix. To construct the secular Hamiltonian matrix we have 
developed rules for evaluating the matrix elements of the h and W operators. Due to the 
generality of the consideration, the developed rules are, of course, applicable for any one- 
and two-body operators. 



2. Variations with respect to the orhitals {(pi} 

The functional differentiation of the energy functional Eq. ljlUp with respect to the one- 
particle functions 0* results in a system of M coupled integro-differential equations 

5i'^\H\^) 

^^J = + ^^^2\<p2) ■■■ + /XiM|0Af), Z = 1, " " " , M. (18) 

By solving these equations for given fixed values of the coefficients Cfi one obtains the 
respective set of the one-particle functions {(pi}. 

The main purpose of this section is to express these equations in an explicit form. To 
fulfil this goal, we first re- write the expectation value of the Hamiltonian, Eq.Q, in a form 
where all the terms depending on orbitals are explicitly visualized and then apply functional 
differentiations. 
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The expectation value of the Hamiltonian can be rewritten in the following form: 

M ^ M 

{^H\^) = J2 Pijh^J + 2 ^ PijkiWijki. (19) 

i,j i,j,k,l 

Here, pij = |\E') are the elements of the reduced one-body density matrix: 

P(^il^i) = ^ y ^2, ■ ■ ■ rAr)^(ri, ra, ■ ■ ■ , fN)df2df3 ■ ■ ■ dr^ 

M 

= E/^^^<^^*(^")<^^(^"'i)' (20) 

id 

and pijki = (\E'|&J&]6fc&i|\E') are the elements of the two-body density matrix: 

p(n, ^2|^i, ^2) = - 1) / ^*(^i> ^2, ^3, ■ ■ ■ rN)'^{fi,f2, rs, ■ ■ ■ , nv)c?^3 ■ ■ ■ c?nv 



M 

= E P^.H</>*(rl)0*(rl)0,(ri)0,(r2). (21) 

i,j,k,l 

The matrix elements of the reduced one- and two-particle densities can be easily obtained: 



Pi'. 



Pij = ^CnC{;n,-l:n,+l;)J ni{nj + 1), (22) 



pun 

ri 



Pijij = ClCfjninj = {h.fi.j). (23) 

n 

Here we present only the diagonal matrix elements of the two-particle density; the off- 
diagonal ones are collected in Appendix EI 

Inspecting Eqs. (j22j) and (j23p we see that the matrix elements of the one- and two-particle 
densities depend only on the expansion coefficients and do not depend on the one-particle 
orbitals {0j} explicitly. In other words, only the one- and two-body integrals /ijj, Wijki 
appearing in Eq. lfT^ have a functional dependence on the one-particle functions 0*. Hence, 
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the functional differentiation has to be apphed only to the integrals hij, Wiju which we treat 
individually according to the general rules of functional differentiation: 

-TT^ = = Wji\(t)k), I ^ J, 25 

Hi Hi 

where we introduce the notation 

Wji = j (t)*{r)W{r-r)(t)i{r)dr (26) 

for the local operators Wji. 

Using Eqs. ()24l25l26|) the variation of the expectation value of the Hamiltonian, Eq. lfTUI) . 
with respect to the 0* takes on the following very compact and appealing form 

^ ' ' ^ = Y.'^Prih + E P^^^i^,,]\ct>,). (27) 

^« 3=1 k,l=l 

Finally, the functional differentiations of the Lagrange energy functional Eq. ()l()|l with respect 
to the one-particle functions result in a system of M coupled integro-differential equations: 

MM M 

Y,Wh + PikjiWkMj) = « = 1, ■ ■ ■ , (28) 

j=i k,i=i j=i 

The main result of the present section is as follows. For any given set of the expansion 
coefficients C^, the best, i.e., energetically most favorable set of one-particle functions {0j} 
used to construct the many-body basis functions (permanents) is determined by solving 
Eq.(|2HI). 



C. Properties of the MCHB(M) equations 

In the formulation of MCHB(M) we assumed that the expansion Eq.® spans a// possible 
configurations of bosons over M one-particle functions. However, the derived equations 
for the optimal expansion coefficients ()12|) and orbitals ()28|) are general and remain valid 
even in the case when the expansion Eq.(jH)) is incomplete and comprises any limited subset 
of configurations. Let us first consider the simplest limiting case, where the expansion Eq. (jS)) 
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contains only a single permanent in which all bosons reside in one and the same orbital 
0i: 

\^) = C^n,)\ni). (29a) 

Of course, rii = N due to the conservation of the total number of particles and C(„j) = 1 
due to the normalization of Then, the expectation value of the Hamiltonian Eq. ()19|l 
takes on the simple form 

(^1^1^) = pn/^ii + ^PiiiiW^iiii (29b) 

and Eq.(j2ni) now reads 

{pnh + PiiiiW^ii} |0i) = yun|0i). (29c) 

The respective non-zero elements of the reduced one- and two-body density matrices 
Eqs. (l22|23j) become trivial 

Pu = {ni) = N, 



Pun = {nl-n^) = N{N -I). (29d) 

Obviously, for contact inter-particle interaction W{f— -P) = XoS{f— r') and putting (j)i = ip 
we easily reproduce the famous GP mean-field energy functional 

Egp = = N[h,, + ^^Ao j ^{f)^d^ (29e) 

and the GP equation for the optimal orbital 

-|^V?^+ V(f) + UN - l)|¥.(f)p| \v) = f^cM, (29f) 

where ficp = P'li/N . 

Next, let us demonstrate that in the more general one-configurational case, when the per- 
manent is given by Eg . (1^ . the many-body MCHB(M) theory boils down to the multi-orbital 
BMP theory 0, 21. |22|. In this case the only permanent contributing to the expansion 
Eq.® with C{n^^n2,- ,nM) = 1 represents a configuration with rii bosons residing in 0i, n2 in 
02, ••• in (pM- Let us first rewrite the general expression for the total energy Eq.()19|) in 
a form where the diagonal Emf and off-diagonal Emb contributions are separated: 

{m\H\^) = Emf + Emb (30a) 
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where 

M ^ M 



""^ 1 1 - % 

Emf = ^[Pii/ijt + ■::PnnWiiii + - ^ {Pijij^ijij + Pijji^ijji)], (30b) 



2'^ 2 
and 

M , M 



-Emb = ^ (Pij/iii + ^ X] PijfciW^ijfcz), (30c) 

where fc,/} runs over all possibilities not included in the respective Emf part. Since 
only one configuration forms the many-body expansion, only the diagonal matrix elements of 



the reduced one- and two-body density matrices given in Eqs. (j22l23|) have non-zero values: 

Pa = {fii) = fii, 
Piiii = {nl - fii) =n^i- Hi, 



Pijij = Pijji = {fiifij) = riiUj. (30d) 
Consequently, Emb = and only the Emf part of the total energy survives and boils down 



to 



M 



M 



(30e) 



i=l 



Analogously, Eq, 

[riih + nj(ri. 



determining the optimal orbitals now reads: 



M 



M 



M 



,)+ ^ ninjWji\(j)j) = ^fiij\(j)j),i = l,--- ,M. 



j=l,j^i 



j=l,j^i 



(31) 



These equations coincide with the multi-orbital mean-field equations, formulated and applied 



in Refs.^ 
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23l|24|. 



We conclude that in the limiting cases where only one permanent forms the many-body 
expansion Eq.(jHJ, the many-body MCHB(M) theory boils down to the respective mean- 
fields. On the other hand, if only one configuration in the many-body expansion is dominant, 
then the multi-orbital mean-field predictions are very close to the many-body ones. Such a 
situation can be realized in real physical systems, for example in deep multi-well traps. 

When the expansion Eq.(jS)) spans all possible configurations of N bosons over M one- 
particle functions, the MCHB(M) equations possess some interesting and useful properties. 
Suppose we solve this system of equations and find the optimal orbitals {(pi} and correspond- 
ing expansion coefficients {C'(ni,n2,n3,---,nM)}- -^y applying an unitary transformation on the 
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0j we can construct a new set of the one-particle functions and find the corresponding new 
set of the expansion coefficients {C{ni,n2,n3,- ,nM)}- This unitary transformation does not 
change (up to a phase factor) the many-body wave-function Eq.®. Consequently, the total 
energy |\E') of the system is invariant with respect to any unitary transformation of the 
one-particle functions. Moreover, we can demonstrate that this is also valid at any iteration 
during the described procedure leading to the self-consistent solution of the MCHB(M) equa- 
tions. The considered invariance of the equations is explained by the fact that the expansion 
Eq.® spans a// possible configurations, i.e., it is complete within the provided subspace of 
one-particle functions. Clearly, for an incomplete expansion this invariance is, in general. 
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lost. For instance, it was shown that the multi-orbital mean-field equations 
which, as discussed above, can be considered as one-permanent MCHB(M) equations, in- 
deed are in general not invariant with respect to unitary transformations of the orbitals 

Having established the equivalence of the MCHB(M) solutions with respect to unitary 
transformations, it is natural to find out which set of orbitals can provide additional physical 
insight. Since the energy is invariant, we consider the reduced one- and two-particle density 
matrices. Having at hand a reduced one-particle density matrix one can diagonalize it: 

M 
i 

The eigenvectors ^f'*^ are referred to as natural orbitals and the eigenvalues pi of this matrix 
are the respective occupation numbers. The natural occupation numbers pi can be considered 
as the average numbers of bosons residing in 0f^*^. The natural orbital analysis of the many- 
body solution is often used to characterize the system. The system is condensed when 
only a single natural orbital has a macroscopic occupation. If several natural orbitals have 
macroscopic occupation numbers then the system is called fragmented js^. Of course, the 
natural orbital analysis can be applied to any state, also to excited states. Comparing the 
natural orbitals and occupation numbers of the ground and excited states one can learn on 
the nature of the excited state. 

Let us examine the properties of the natural orbitals in the MCHB(M) theory. On the one 
hand, we know that each eigenvector of the one-body density matrix p(r|F) can be expressed 
as a linear combination of the MCHB(M) orbitals. On the other hand, we have seen that 
due to the completeness of the many-body expansion, the MCHB(M) solution is invariant 
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with respect to any unitary transformation applied. We conclude that the natural orbitals 
themselves constitute a MCHB(M) solution as well. This finding gives us the freedom to 
characterize the MCHB(M) one-particle functions by natural orbitals. We shall do so in the 
following, unless explicitly mentioned. 

Using the multi-configurational expansion Eq.Q as an ansatz for the true many-body 
wave function of a trapped bosonic system and applying the variational principle we derived 
general equations which allow one to find the best many-body basis functions and the cor- 
responding expansion coefficients self-consistently. We have shown that these equations are 
also applicable in any desired active sub-space, including the limiting case of one-permanent. 
In the latter case the derived equations boil down to the one-orbital Gross-Pitaevskii and 
multi-orbital mean-fields equations. At this point it is very important to stress that all 
derived equations and conclusions are general, i.e., they are valid for any geometry of the 
external trap potential, for any physical shape of the inter-particle interaction and in any 
dimension. 



III. ILLUSTRATIVE NUMERICAL EXAMPLES AND APPLICATIONS 



A. Preliminaries and implementation of MCHB(2) 

In the present section we implement the developed MCHB(M) formalism for systems 
of cold bosonic atoms. We consider the simplest case of M=2, i.e., MCIIB(2) theory and 
apply it to the ground and excited states of the bosonic systems trapped in symmetric 
and asymmetric one- and two- dimensional double-well potentials. While the ground state 
of the one-dimensional bosonic gas trapped in symmetric double-well potentials has been 
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35| and references therein, the 



intensively studied in the literature see Refs. 
many-body properties of bosonic systems at higher dimensions in symmetric and especially 
asymmetric double-well traps Q] are open theoretical questions. 

The configurational space of the MCHB(2) theory constitutes all symmetrized permuta- 
tions of N bosons over two orbitals 0i and 02- The trial many-body function Eq.Q in this 
case is spanned by N+1 configurations: 

N 

1^) = XI ^ini,n2)\ni,n2) (33) 

711=0 
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The condition rii + n2 = N implies that only one occupation number rii can be used to 
enumerate the configurations. 

The configurational space of MCHB(2) theory also implies that the non-zero Hamiltonian 
matrix elements can be only between the three generic permanents 

\Po) = 1^1,^2), |Pi) = |ni + l,n2-l), IP2) = |ni + 2,n2-2). (34) 

Using these permanents and the general rules for evaluating one-body and two-body matrix 
elements, Eqs. fll4ll5ll6ll7j) . the expectation value of the Hamiltonian (fT^ becomes 

{%H\^) = {ni)hn + {n2)h22 
^{r^i^Mw.n, + - ^^^^ 2222 + {nin2)iW,2u + Wu2i) 
+ Puhu + P2ih2i + P2iiiW^2iii + PnuWnu 

+ P222lW^2221 + Pl222W^1222 + — M^2211 + — W^1122- (35) 

The coupled system of integro-differential equations needed for determination of the 
optimal one-particle functions reads: 

+ {{nl) - {ni))Wn + (^1^2)1^22 + P2iiiW^2i + PuuWu] |0i) + 

|(nin2)W^21 + PnuWii + pi2h + Pl222W^22 + Pll22W^12| |02) = + Pl2|02) 

|(n2)/l + {{flD - (n2))#22 + {nih2)Wii + P222lW^21 + Pl222W^12} \<p2) + 

|(fiin2)#i2 + P222lW^22 +P2l/i + P211lW^ll + P221lW'2l} \<Pl) = fJ'2l\4>l) + P22\4>2) ■ 

(36) 

As an illustrative example we consider N identical spinless bosons, interacting via contact 
potential iy(r*— F) = Xo6{f— -P), where Aq is proportional to the s-wave scattering length. 
In this case the two-body integrals, see Eq.®, appearing in Eq. (j35|) simplify to 

and instead of the local integral operators appearing in Eq. ()36|) we obtain the functions 

This considerably simplifies the implementation of the MCHB theory as the system of 
integro-differential equations (|36p boils down to a system of non-linear coupled differen- 
tial equations. Furthermore, in this paper we work in = 1 units, where m is the mass 
of boson and L is the length scale. 

20 



We recall that for MCHB(2) theory the number of particles, the inter-particle interaction 
strength Aq and the external trap potential have to be specified. The expansion coefficients 
{C(ni,n2)} orbitals 0i and 02 are treated as variational parameters. The optimal values of 
these parameters are determined self- consistently by diagonalizing the secular Hamiltonian 
matrix to find the expansion coefficients, Eq.(|12j). and by solving the coupled system of 
non-eigenvalue non- linear differential equations Eq.()36|) to get the one-particle functions. 

The self-consistent procedure of finding the optimal one-particle functions and corre- 
sponding expansion coefficients can be implemented as follows. We start from some initial 
guess for the one-particle functions (pi, 02 obtained, say, as the two eigenfunctions lowest in 
energy of the bare Hamiltonian 

1, 



(37) 



where i = 1,2 and V{r) is the corresponding trap potential. These one-particle functions are 
used to construct the permanents in the expansion Eq. ()33j) and to evaluate the Hamiltonian 
matrix elements. By diagonalizing this matrix, we solve the corresponding secular equation 
Eq.()12j) and get N+1 orthonormal eigenfunctions and the corresponding eigenvalues. The 
lowest-energy solution corresponds to the ground state of the problem, while other solutions 
may be used to attack excited states. The eigenvector of interest contains the set of expansion 
coefficients {C(„^^„2)}. To implement the second part of the variational principle, we use these 
expansion coefficients to evaluate the elements of the reduced one- and two-body densities 
given in Eqs. ()22l23|) and in AppendixEl required by the coupled system (jHUj) of the MCHB(2) 
equations. By solving this system of equations we find a new set of the one-particle functions 
01,02- This self-consistent scheme is iterated until convergence is achieved. 

One goal here is to demonstrate that when the bosonic system is treated at the many-body 
level and the many-body basis set is incomplete, then the choice of one-particle functions 
used to construct this many-body basis set has a great impact on the results and predic- 
tions obtained. To realize this in practice we consider two sets of one-particle functions and 
compare the many-body predictions obtained within each set. As the first set we use the 
lowest-energy eigenfunctions of the bare Hamiltonian Eg. (13711 . We stress that such a fixed 
choice of the one-particle functions is often used [3, [28^ in many-body treatments on 
bosonic systems. The MCHB(2) solution is the second set of orbitals which is, of course, the 
best possible choice because these orbitals have been obtained in the framework of the full 



21 



variational principle. Prom now on we use BH and MCHB superscripts to distinguish the 
results obtained within bare Hamiltonian and MCHB(2) one-particle function sets, respec- 
tively. For example, we denote the corresponding one-particle functions as 0f ^, (/>f ^ and 
^MCHB ^ ^MCHB ^ In GUI iterative MCHB(2) scheme we use 0f^,0f'^ as the initial guess. 
Therefore, by comparing the results obtained at the first and last iterations we immediately 
observe the impact of self-consistency. 

Let us now elaborate on the choice of trapping potentials. In our work we have chosen 
the following form of symmetric and asymmetric double- well potentials. Let us imagine two 
separate trap potentials, for simplicity in ID, described by Vi{x) and V2{x). We construct 
a "superposition" of these traps in such a way that, on the one hand, the profiles of each 
potential well in the resulting double-well potential would be as much as possible close to 
the original potentials. On the other hand, it is also desirable that the inter-trap separation, 
degree of asymmetry and barrier height can be easily manipulated. What is also important 
is that this double-well potential would have a simple (differentiable) analytic form and 
permit 2D and 3D generalizations. Let us diagonalize the matrix 



Vi{x + xo) + B C 

C V2{x-xo) 



(38) 



The lowest eigenvalue of this matrix is a function which fulfils almost all the above men- 
tioned requirements. The original traps are connected to each other and the parameter C 
is responsible for the smoothness of this connection. The minima of the wells are located 
approximatively at ±xo and the Vi{x) trap is biased by B with respect to the V2{x) one. 
However, to gain an additional control on the barrier height we add a smooth barrier function 
centered at the origin: 

where D defines the width of this additional barrier function and A can be used to vary its 
height. The resulting trap potential takes on an analytic form: 



y{^) = -[Vi{x + xo) + B-^jAC'' + {Vi{x + xo) + B- V2{x - xo)y + V2{x-xo)]+Vb{A, D). 



We recall that in this paper we work in = 1 units, where m is the mass of boson and L 
is the length scale. 
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B. Ground state in ID symmetric double-well trap 



Let us first apply MCHB(2) theory to study the ground state of repulsive bosons trapped 
in a symmetric double- well potential. This is a popular problem, widely discussed in the 



literature, see Refs.bO 
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35 1 and references therein. We consider a system of N=1000 

□ 

bosons which is of the order of the particle number taken in recent experiments j9||. This 

system is trapped in the symmetric double-well potential: 

1 



Vsyn^M = 0.5x2 _ _^25 + 64a;2 + 8.19531 + Vb{A = 4,D = 0.75). (40) 

This potential plotted in the left panel of Fig^ by a thick solid(black) line is obtained by 
diagonalizing Eq.§^ with Vi{x) = 0.5{x + 4.0)^, V2{x) = 0.5(x - 4.0)^, C = 2.5, B = 0.0 
and by adding a barrier function Vb{A = A,D = 0.75) (see Eg. ()39|l ). A constant shift has 
been applied to put the minimum of the potential energy to zero. 

In the right panel of FigC] we plot the ground state energies per particle of this system 
as a function of the inter-particle interaction strength Aq. The many-body results obtained 
within fixed bare Hamiltonian (BH) one-particle functions are plotted by a solid (black) 
line with open circles. The solid (red) line with filled circles represents the energies per 
particle obtained self-consistently within the framework of the MCHB(2) theory. Both curves 
coincide only at the limit of non-interacting particles, the optimal one-particle functions in 
this case are the bare Hamiltonian functions. Increasing Aq, the MCHB(2) energy curve 
gradually (exponentially) deviates from the fixed-orbital one. To illustrate this we plot 
in the inset the difference between these energies per particle A = E^^ — E^'^^^^ as a 
function of Aq (notice the logarithmic scale on both axes). It is clearly seen that the many- 
body treatment with fixed one-particle functions provides an adequate description only up 
to Ao ~ 10^^. For stronger inter-particle interactions the self-consistency becomes more and 
more relevant, and leads to lower energies. We conclude that the choice of the one particle 
functions in truncated CI expansions is very crucial for the appropriate description of the 
energetics of interacting bosonic systems. 

Next, we construct the reduced one-particle density matrix Eq. ()22|l and find the corre- 
sponding natural orbitals and natural occupation numbers Eq. ()32|l . The many-body ansatz 
used in the MCIIB(2) theory implies that there are only two natural orbitals with respec- 
tive occupations numbers pi and p2- The conservation of the total number of particles 
Pi + P2 = N allows us to consider only one occupation number. 
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In FigI21 we plot the occupation number p2 of the second natural orbital as a function 
of the inter-particle interaction strength Aq. The solid (red) line with filled circles and the 
solid (black) line with open circles represent the self-consistent MCHB(2) and fixed-orbital 
BH results, respectively. Up to Aq ~ 10^'^ both methods predict that the value p2 gradually 
increases with Aq. In other words, according to both many-body treatments fragmentation 
takes place when Aq is increased and at Aq ~ 10"'^ approximately 200 particles out of 

= 1000 are fragmented. However, for Aq > 10^^ the predictions obtained within the 
self-consistent MCHB(2) and the fixed-orbital many-body theory start to deviate from each 
other and eventually drastically contradict each other. The many-body results obtained with 
fixed bare Hamiltonian functions show that the fragmented fraction p2 increases further 
with Aq until it saturates to some constant value, p2 ~ 366. In contrast, the MCIIB(2) 
theory predicts that at some inter-particle interaction strength the p2 fragmented fraction 
approaches its maximum value and then gradually decreases. Finally, we would like to 
mention that at much larger values of Aq another physical phenomenon starts to take place, 
fermionization but the region of inter-particle interaction strengths studied here is 

far below this limit. 

Now we elaborate on the physics behind these many-body predictions. The ground state 
fragmentation phenomenon studied here appears due to the double- well topology of the trap 
potential and disappears at zero barrier height. With increasing inter-particle interaction 
the respective chemical potential(s) of the trapped repulsive bosons increase as well and this 
can be viewed as an effective reduction of the barrier height. The trap potential used in 
the present study has a barrier of finite height and, hence, from some critical interaction 
strength on the bosons are energetically above the barrier and do not "see" it. We may thus 
conclude that the self-consistent MCIIB(2) theory predicts a physically relevant behavior 
of the fragmentation as a function of inter-particle interaction strength in contrast to that 
predicted by the fixed bare Hamiltonian functions. 

Now we investigate the fragmentation phenomenon in the symmetric double-well potential 
as a function of the barrier height. We ask the question at which barrier height the system 
of = 1000 bosons at fixed inter-particle interaction strength of Aq = 0.01 becomes, say, 
25% fragmented. We consider the same symmetric double- well trap potential as in Eq. ()40|) 
and rump the barrier up by varying the parameter A of Vb(y4, D = 0.75) (see Eq.(jnH)). For 
every A a constant shift is applied to put the minimum of the respective potential energy 
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to zero, hence Vsymm{x) at x = determines the value of the barrier height. In Fig|21we 
plot the occupation number p2 of the second natural orbital as a function of the barrier 
height. The solid (red) line with filled circles and the solid (black) line with open circles 
represent the MCHB(2) and fixed-orbital BH results, respectively. In this figure it is clearly 
seen that 25% fragmentation ratio, i.e., p2 = 250 out of = 1000 bosons, is obtained in 
the framework of fixed-orbital many-body theory at a barrier height Vsymm{0) ~ 6.75, while 
the self-consistent MCHB(2) gives such a fragmentation ratio when the barrier height is 
Ks/mm(0) ~ 10.3. Again, the BH predictions considerably overestimate the fragmentation. 
In reality the fragmentation develops slower with increasing barrier height than predicted by 
the fixed-orbital BH many-body method. This observation is also of a relevance for multi- 
well systems, including optical lattices. We stress that the difference between predictions 
of both theories has a non-trivial origin, and one curve can not be obtain from the other 
by a simple procedure, e.g. shift. To illustrate this, we plot in FigOlby dashed (blue) line 
the difference between corresponding occupation numbers A = — p^^^^^ as a function 
of the barrier height. In this figure we see that this difference is substantial at any finite 
barrier heights and becomes less pronounced only in the limit of very large barrier heights. 

In these investigations of the ground state of a bosonic system trapped in symmetric 
double-wells we have seen that the predictions obtained within the framework of the fully 
self-consistent MCHB(2) and within fixed-orbital many-body theories utilizing the same 
active space are quantitatively and some times even qualitatively different. By construction, 
the self-consistent description is more precise than the fixed-orbital one. The fixed-orbital 
many-body theory can, in principle, reproduce the self-consistent MCHB(2) results if more 
BH orbitals are included, i.e., if the active space is enlarged, resulting in a substantial 
increase of the computational effort which can be in practice beyond reach. 



C. Ground state in ID asymmetric double- well trap 
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Despite considerable progress in the experimental studies on double-well traps jS 
a realization of a double-well potential with perfect symmetry is not straightforward. In 
contrast, bosonic systems trapped in a perfect double-well pot ential is the most popular 
theoretical problem addressed in the literature 
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on bosonic systems in asymmetric traps remain scarce 



I, while theoretical studies 
To elaborate on this very 
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complicated problem, we address here the ground state of bosonic systems trapped in one- 
dimensional asymmetric double-well trap potentials. Again, the goal is to demonstrate that 
self-consistent many-body methods remain of importance. 

To construct the asymmetric double- well trap we locate two equal harmonic traps at ±xo 
and displace the left trap upwards with a bias B. Clearly, if the wells are well separated, then 
the bare eigenfunctions of this trap lowest in energy are predominantly localized either in the 
left or right wells and keep the shapes of the pure harmonic functions. For a comparatively 
small bias B, the three eigenstates of the double- well potential lowest in energy are ordered as 
depicted in the left panel of FigI3 Such an asymmetric double- well potential can be obtained 
by diagonalizing Eq.^ with Vi{x) = 0.5{x + A.0f, V2{x) = 0.5(x-4.0)^ C = 2.5, B = 0.5 
and by adding a barrier function Vi,{A = A,D = 0.75). A constant shift is also applied to 
put the global minimum of the potential energy to zero. The analytical function of this 
potential reads 

K,^.„(x) = 0.5x^ - iv25.25 + 8x + 64.= + H(A = 4,D = 0,75) + 8.4423. (41) 

We consider N=1000 bosons trapped in the asymmetric double-well potential of Eq.()4ip. 
Let us see how the ground state of this system develops with increasing inter-particle inter- 
action strength Aq. The physical picture of this development, also supported by mean- field 
studies is as follows. At zero interaction all bosons are localized in the deeper 

(right) well. The bosons continue to stay localized in this well up to some critical inter- 
action strength Ag**. From this Ag^ on the tunnelling of bosons into the left well becomes 
possible and bosons start to occupy the left well. In other words, there are two regimes of Ag, 
in the first regime Aq < Ag*" the ground state properties depend mainly on the geometry of 
the deeper well, while in the second regime Ag > Ag** they depend on the global geometry of 
the asymmetric double- well potential. These observations are supported by our MCHB(2) 
calculations presented in the right panel of Fig0] The MCHB(2) theory gives Ag'' = 0.00136 
for the transition between the two regimes. 

In the right lower panel of FigHJone can see that in the first regime (Ag = 0.00135 < Ag'') 
both the orbitals and the density are indeed localized in the right well. The natural analysis 
tells that pi = 999.98 bosons are condensed on the first natural orbital (red) and the fraction 
of p2 = 0.02 bosons is depleted to the second natural orbital (blue). Since both orbitals are 
localized in the same well, we can say that the origin of the depletion is on-site excitations. 
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The MCHB natural orbitals and the density for the second regime (Aq = 0.01 > Ag'') are 
plotted in the upper right panel of FigE] and as it was expected in this case both wells are 
populated. This ground state is almost 10% fragmented, because pi = 906.06 bosons reside 
in the first and p2 = 93.94 bosons in the second natural orbital. 

Let us see whether it is possible at all to obtain qualitatively the same results by using 
the fixed-orbital many-body method. The bare Hamiltonian functions of the asymmetric 
double- well potential of Eg. 1)411) are depicted schematically in the left panel of Fig^ If 
the first and second orbitals are used to construct the permanents, then at any non-zero 
interaction strength the bosons are spread over both wells. Consequently, with such a choice 
of one-particle functions the first regime can not be described. Instead, one can try to use 
the first and third eigenfunctions of the bare Hamiltonian to construct the many-body basis 
set. In this case, however, it is impossible to address the second regime. To overcome this 
difficulty one can use all three orbitals simultaneously. However the active space, i.e., the 
number of many-body basis functions used in this case is much larger then the MCHB (2) 

/ 3 + 1000 -l\ 

ones. For = 1000 we would need = 501501 configurations instead of 

y 1000 J 

1001! Still, the self-consistency is not used and the quality of these fixed-orbital results has 
to be investigated. 

D. Distributions of the expansion coefficients 

So far, to study bosonic systems trapped in the symmetric and asymmetric double-well 
potentials we have considered and analyzed the MCHB energies and orbitals. We recall that 
the MCHB solution is given by the optimal sets of one-particle functions and of the respective 
expansion coefficients obtained self-consistently. In this section we analyze in some detail the 
MCHB coefficients {C(ni,--- ,nM)} appearing in the expansion Eq.® and exploit the freedom 
of unitary transformation as put forward in Sec lII CI 

Generally, for any orthogonal many-body basis set the square of the expansion coef- 
ficient C^*jj^ ... „^j^C(^ni,--- ,nM) defines the probability to find the many-body solution in the 
configuration described by the respective many-body basis function |ni, n2, ns, ■ ■ ■ ^um)- 
In other words, the squares of the expansion coefficients can be considered as a multi- 
dimensional discrete probability distribution in the discrete sample space spanned by the 
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many-body basis functions. In the MCHB(2) theory, to count all many-body basis functions 
\ni, ^2) = \ni, N — rii) one needs only one independent parameter rti {n2 = N — rii). There- 
fore, the squares of expansion coefficients can be viewed as a probability function P{n) of a 
discrete distribution defined over the independent parameter n = 0,1,2, ■■ ■ ,N: 

Pin) = Cl^j,_^f^n,N-n). (42) 

We can use the mean value v and variance as measures of the distribution: 

N N 

ni=0 rai=0 
N N 

< = E ^(^On? - ( Yl nni)n,f ^ {hi) - {h,)\ (43) 

?ii=0 ni=0 

Here, we use rii as the independent parameter. If the occupation number of the second 
orbital ^2 is used instead then = N — Interestingly, the {hi) and {hf) have already 
appeared in the evaluation of the diagonal elements of the reduced one-particle Eq. and 
two-particle Eq.(j2Sl) density matrices. 

As mentioned above in Sec lII CI due to the invariance of the MCHB solution with respect 
to unitary transformations, there are infinitely many possible choices of the MCHB orbitals 
which give the same energy. It is also clear that the distribution of the expansion coefficients 
depends on the particular choice of the one-particle functions used to construct the many- 
body basis set. Therefore, there are infinitely many possible distributions of the expansion 
coefficients as well. However, since the mean values and variances of the distributions are 
different, we use these quantities as the main criteria to compare energetically equivalent 
distributions. The main aim now is to find the distributions characterized by the minimal 
width. 

Let us consider two sets of the MCHB(2) orbitals connected by a unitary (orthogonal) 
transformation 

cos 9 sm.6 \ I (hi \ , , 

' (^^) 

— sin 9 cos 9 j \ 4>2 I 

where 9 is the rotation angle. Clearly, the creation and annihilation operators corresponding 
to each set of the orbitals are coupled via U as well: 













[hi 









\ I cos 9 sin 9 

Li] 
h 



bn^ I \ — sin 9 cos 9 




(45) 



28 



Since all possible real MCHB(2) orbitals can be obtained from the initial (0i,02) ones by 
changing the angle 9, the respective distributions of the expansion coefficients as well as 
their mean values and variances depend also on this angle. Here we are interested in the 
variance: 

~< = (^?) - (^i)' = (h\hih\h,) - {b\hf = a\e). (46) 
After straightforward algebra one finds 

+ sin^ 6 cos^ 9[A{hih2) - 2{hi) (^2) + + P1122 + P2211 - (pi2 + P2if] 
+ 2 sin 6^ cos^ 9[pi2 + P1112 + P2111 - (Pi2 + P21) (^1)] 

+ 2sin^6'cos6'[p2i +P1222 + P2221 - (Pi2 + P2i)(^2)], (47) 

where the involved diagonal and off-diagonal elements of the reduced one- and two-body 
density matrices are given in Eqs. ()22l23|) and in Appendix ^ The extrema of this function 
are obtained in the ordinary way, by solving 

^a^e) = 0. (48) 

The procedure of finding the distribution with minimal variance is implemented as follows. 
Having at hand a MCHB(2) solution, i.e., the orbitals (0i,02) and a set of the respective 
expansion coefficients {C(^ni,n2)}: we recompute all required elements of the reduced one- 
and two-body density matrices, appearing in Eq. ()47|) and explicitly reconstruct the variance 
(T^(0). The angle 6min at which this function has a minimum is obtained numerically by 
solving Eq. ()48|l . The unitary transformation Eq. ()44|1 with this angle gives a new set of 
MCHB(2) orbitals (0i,02)- The distribution of expansion coefficients computed on these 
orbitals has the minimal variance amm- Obviously, we can call such a distribution the 
minimal distribution. 

The natural orbitals being the MCHB solutions can be used to shed light on the physical 
nature (depletion or fragmentation) of the ground and excited states. Let us see how the 
expansion coefficients distributions of the many-body basis sets constructed by using the 
natural orbitals look like. In the left lower panel of Fig El we plot the expansion coefficients 
obtained for the ground state of N=1000 bosons with Aq = 0.01 trapped in the symmetric 
double- well potential Eq.lflUI). This system has been discussed above in Sec lIIIBl The 

29 



corresponding ground state natural orbitals are very similar to those plotted in the right 
lower panel of FiglHl Due to the symmetry of the trap potential, the ground many-body 
state is, of course, of gerade symmetry, while the natural orbitals are of gerade and ungerade 
symmetries. Therefore, for the ground state the non-zero expansion coefficients can appear 
only due to the contributions of configurations of gerade symmetry. Indeed, in the left lower 
panel of Fig0 one can see that the configurations with even number of bosons residing in 
the ungerade orbital contribute, while those with odd numbers i.e., |1000 — 1, 1), |1000 — 
3, 3), 1 1000 — 5,5) etc., give zero contributions. In this figure it is clearly seen that the 
main contributions come from the first configurations, where almost all bosons reside in the 
first orbital, while the configurations with large populations of the second orbital do not 
contribute. This observation is supported by the statistical description of this distribution 
in terms of its mean values (see Eq. ()43|) ). The mean statistical values of this distribution 
are i>i = (ni) = 994.78 and z/2 = (^2) = 5.22. As one would expect these mean values are 
identical to the natural occupation numbers of the respective natural orbitals. 

In the right lower panel of FigEl we plot the expansion coefficients obtained for the 
system of N=1000 bosons with Aq = 0.01 trapped in the asymmetric double-well potential 
Eq. (PT|) discussed in Sec lIII Cl The corresponding natural orbitals presented in the right 
upper panel of FigE] do not possess any symmetry. Consequently, all the many-body basis 
functions constructed by using these natural orbitals can give non-zero contributions in the 
many-body expansion. Indeed, we can see this in the right lower panel of FiglHl The mean 
statistical values of this asymmetric distribution ui = (fii) = 906.06 and 1^2 = (^2) = 93.94 
are, of course, identical to the respective natural orbital occupation numbers. 

Let us see how the variance cr^^ = (ni) — (ni)^ characterizes the distributions of the 
expansion coefficients. Inspecting the distributions obtained with natural orbitals we see 
in the lower panels of Fig|Sl that for the symmetric double- well, where only several many- 
body basis function have non-zero expansion coefficients, the variance is quite small a no = 
7.654. While in the asymmetric case, where almost all expansion coefficients have non-zero 
contributions, the width of the distribution is much larger a^o = 121.731. We conclude that 
the value of the statistical variance cx^, characterizing the width of the distribution of the 
expansion coefficients, indeed provides an adequate estimation on the number of contributing 
configurations. 

Now let us see how the minimal distributions, i.e., the distributions with minimal width 
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look like. We consider the same examples of the symmetric and asymmetric double- wells as 
above, and find their minimal distributions by using the developed algorithm (see Eq. ljlHjl ). 
The minimal distributions obtained for the symmetric and asymmetric double-wells are 
presented in the left and right upper panels of FigEl respectively. The respective variances 
are also depicted. The minimal distribution obtained for the symmetric double-well trap 
has a maximum located exactly at the |500,500). This is a symmetric distribution because 
the pairs of basis vectors around maximum |500 — i, 500 + i) and |500 + i, 500 — i) contribute 
with identical coefficients. Interestingly, the minimal distribution is smooth, in contrast to 
that obtained within the natural orbitals, where the neighboring MCHB(2) coefficients are 
of different sign, see lower left panel of FigO The width of the minimal distribution is, of 
course, smaller than that for the natural orbital, = 3.313 in comparison to a^o = 7.654. 

Comparing the right upper and lower panels of FigjSlwe see that in the presented exam- 
ple of asymmetric double-well trap the minimal distribution differs substantially from that 
obtained within the natural orbitals. For the asymmetric case, the width of the minimal 
distribution amin = 0.753 is about 2.5 orders of magnitude smaller than that for the natu- 
ral orbital crjvo = 121.731. This means that the minimal distribution is formed by several 
configurations in contrast to the broad distribution obtained with natural orbitals where 
all configurations contribute. The maximum of the minimal distribution is located around 
|600,400). The distribution is smooth and, of course, not symmetric. 

Irrespective of the symmetry of the trap potential used, the width of the minimal distri- 
bution can be much smaller than that obtained with the natural orbitals. The profile of the 
minimal distribution is smooth, while the distribution of the expansion coefficients obtained 
within other orbital sets is not necessarily so. These observations lead to several important 
consequences. First, the minimal distribution of the expansion coefficients can be approxi- 
mated by a smooth continuous function. Looking at the pictures in Fig^we approximate 
the probability function Eq. ()42|l of the minimal distribution by a Gaussian function: 



The parameters of this function are obtained straightforwardly. We take ^ = rii as an 
independent variable. The averages Eq. ()43|1 of the minimal distribution define the location 
■Co = (^i) and width a = (Tmin of the Gaussian. In the upper panels of FigEl we plot 
the Gaussian distribution functions approximating the minimal distributions of the studied 
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symmetric and asymmetric systems by black solid lines. We stress that in Fig|31 we plot 
the distributions of the expansion coefficients, i.e., a/ P{ni). In this figure we see that the 
Gaussian distributions match the numerical results very well. Moreover, once a MCHB 
calculation has been performed, this continuous Gaussian approximation does not require 
any fitting parameters. 

We mention that continuous approximations to the discrete distributions of the expansion 
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coefficients have already been addressed in the literature for symmetric trap potentials 
I40I . In these studies smoothness of the real CI coefficients and thus of the continuum 
distribution approximation was used as a basic assumption. In contrast, here we have 
demonstrated numerically that the profile of the distribution of the expansion coefficients 
obtained within MCHB many-body method can indeed be smooth. Moreover, the smooth 
profiles of the discrete distribution of the expansion coefficients are observed in examples of 
symmetric as well as of asymmetric trap potentials. However, it is very important to stress 
that smoothness is achieved only for the minimal distributions, i.e., for a very specific choice 
of the orbitals (see Eq. ljTfj) ). 

Finally, we remark that the existence of smooth continuous functions approximating the 
discrete distribution of the expansion coefficients makes the developments of self-consistent 
many-body methods very promising for attacking many-particle bosonic systems within huge 
configurational spaces. 



E. Excited states in ID symmetric double-well trap 

In the preceding sections we considered the ground state of a trapped bosonic system 
and addressed condensation and fragmentation as properties of the ground state. The main 
subject of this part of our work is to touch upon properties of excited states of trapped 
bosons. The studies on excited states are of great interest Q Q, Q Q, Q because of 
their relevance for depletion and stability of condensates, for time-dependent and finite- 
temperature effects, for formations of solitons and soliton trains js^, as well as for other 
interesting phenomena. 

Let us assume that we have obtained the self-consistent MCHB orbitals for the ground 
state of the bosonic system. Since in the MCHB scheme the diagonalization of the secular 
matrix was employed we also have the energies and many-body wave-functions of the excited 
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states. However, the excited states computed in this way are not self-consistent. Here we 
address the question whether the self-consistent orbitals obtained for the ground state can 
also be used to provide an adequate description of the excited states or whether every excited 
state has to be treated individually. 

To answer this question we consider a system of N=1000 bosons with Aq = 0.01 trapped 
in a symmetric double-well potential. In this example we use the trap: 

Vsy^m{x) = 0.5x^ - iV25 + 64x2 + Vb{A = 25,D = 0.75) + 8.19523, 

obtained by diagonalizing Eq.^ with Vi{x) = 0.5(x + 4.0)2, V2ix) = 0.5(x-4.0)^ C = 2.5, 
B = 0.0 and by adding a barrier function Vh{A = 25, D = 0.75). To put the minimum of 
the potential energy to zero we also use a constant shift. 

First, we apply the MCHB(2) approach to obtain the self-consistent energy and orbitals 
of the ground state. This state is essentially condensed because the occupation numbers of 
the corresponding reduced one-particle density matrix are pi = 994.78 and p2 = 5.22. In 
other words, 994.78 particles are condensed in the first natural orbital and 5.22 are excited 
out to the second orbital. The density and respective natural orbitals are plotted in the 
lower right panel of FiglHl We recall that the natural orbitals are solutions of the MCHB. 

Having at hand the self-consistent ground state orbitals we diagonalize the full secular 
matrix and get the energies of the excited states. In the left panel of FiglHl we depict the 
energy level of the first excited state. We connect both levels by a vertical solid line with 
arrow to stress that this first excited state is obtained by using the MCHB orbitals of the 
ground state. The natural orbital analysis applied reveals that the occupations numbers 
of this first excited state are pi = 985.20, p2 = 14.80 and the natural orbitals of this 
excited state are almost identical to the ground state ones. Comparing the natural orbital 
occupation numbers of the ground and this first excited state we conclude that this excited 
state can be considered as a further microscopic excitation of a small number of particles 
out of the condensate, i.e., a further depletion of the condensate. 

Let us now see what happens when the first excited state is treated self-consistently. To 
realize this we employ the developed MCHB(2) method to the excited states as follows. 
We recall that in our iterative MCHB scheme the many-body expansion coefficients corre- 
sponding to the ground state are obtained as components of the first, i.e., lowest-energy 
eigenvector of the secular matrix. To attack the first excited state we use the components of 
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the second eigenvector during all iterations. The self-consistent results obtained by applying 
this procedure to the first excited state are also depicted in FiglHland marked as "SC Excited 
State". The natural orbitals and the density corresponding to this state are shown in the 
upper right panel and its energy level is depicted in the left panel. 

From Figiniit is clearly seen that self-consistency can have an enormous impact on prop- 
erties of an excited state. The energy of the first excited state obtained self-consistently is 
much lower than that obtained by using the ground state orbitals. The self-consistent first 
excited state is almost 30% fragmented. The respective natural orbital occupation numbers 
are pi = 734.84 and p2 = 265.16, contrasting pi = 985.20 and p2 = 14.80 obtained above 
for the non-self-consistent excited state. By comparing the upper and lower right panels of 
Figiniwe see that the shapes of the first natural orbitals (red) in both calculations are almost 
identical, while those of the second natural orbitals (blue) differ drastically from each other. 
Moreover, it is easily seen in the right upper panel of Fig IHl that a simple linear combination 
of the natural orbitals of the first self-consistent excited state gives almost pure left and right 
localized functions, which are, of course, solutions of the MCHB as well. For the ground 
state orbitals depicted in the right lower panel of FiglHlsuch a localized presentation can not 
be obtained. From this analysis we conclude that the first self-consistent excited state and 
the ground state are qualitatively different and exhibit a different "topology". 

In this example the ground state of the system is condensed and the first excited state is 
fragmented. Fragmentation is shown to be much more favorable energetically than a further 
depletion of the condensate. One goal of this study is to demonstrate that self-consistency 
can be very important for an adequate description of excited states. Indeed, we have shown 
that without self-consistency the lowest-energy excited state describes a depletion of the 
condensate. The inclusion of self-consistency significantly lowers the energy of the first 
excited state and drastically changes its character. Instead of depletion of the condensate 
it describes its fragmentation. We stress that excited states of different topology also exist 
in many other trapped bosonic systems. The question whether they are low-lying or highly 
excited states depends on trap geometries, number of particles and strength of inter-particle 
interaction. 
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F. Ground state in 2D symmetric double-well trap 

In the present section we investigate the relevance of self-consistency for many-body 
studies on trapped bosonic systems in higher dimensions. For illustration, here we investigate 
a repulsive bosonic system trapped in the two-dimensional symmetric double-well potential 

Vsy^U^, y) = 0.5x^ + 0.5y^ - ^V25 + 64x2 + 8.19531 + Vb{A = 8,D = 0.75). (50) 

This potential is obtained according to Eq. ljHHj) as a superposition of two pure harmonic 
2D potentials Vi{x,y) = 0.5[(x + 4.0)^ + y^] and V2{x,y) = 0.5[(x - 4.0)^ + y"^] with 
C = 2.5, B = 0.0 and by adding the two-dimensional barrier function Vb{A, D) = 
exp [—{x + yY/ {20"^)]. Here, we have also applied a constant shift to put the min- 
imum of the potential energy to zero. 

The ground state of N=1000 identical bosons at Aq = 0.01 trapped in this double-well 
trap has been investigated within the framework of the fixed-orbital and self-consistent 
MCHB(2) many-body methods. The two eigenf unctions lowest in energy of the respective 
two dimensional bare Hamiltonian Eg. ()37|1 have been used to construct the fixed-orbital 
many-body basis set. As in the one-dimensional case, we use these functions as an initial 
guess for solving MCHB(2) equations. 

The geometry of the double-well trap used implies that the ground state density is made 
of two parts each localized in one well. Due to the perfect two-fold symmetry of the trap 
potential it suffices to consider only one of them without loss of information. In Fig|3 for 
convenience of comparison, we depict the part of the self-consistent MCHB(2) ground state 
density localized in the left well together with the part of the total density obtained within 
the usual fixed-orbital many-body method localized in the right well. From this figure it 
is clearly seen that the densities obtained are very different. To better account for the 
repulsion between the bosons, the self-consistent density is more delocalized within the well 
than the fixed orbital one. Of course, the MCHB(2) solution has a lower energy than the 
fixed orbital one. The natural orbital analysis applied shows that in the MCIIB(2) case 
pi = 750 bosons reside in one orbital and p2 = 250 bosons in the other orbital, while for 
the many-body solution obtained with fixed bare Hamiltonian orbitals these occupations are 
Pi = 634 and p2 = 366. Both many-body methods give the same qualitative prediction on 
the nature of the ground state - this state is fragmented, however, the predicted details of 
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the fragmentation are very different. 

In this investigation we have demonstrated that analogously to the ID case, self- 
consistency is of great relevance for the many-body description of bosonic systems in higher 
dimensions. 

IV. IMPLEMENTATION OF THE MCHB(M) FOR TWO BOSONS IN A TRAP 

The technical realization of the developed MCHB method for the general case of M 
orbitals and N bosons requires considerable methodological and algorithmical efforts. In 
this section we perform the first step in this direction and implement the MCHB(M) theory 
for two interacting bosons. 

We consider a system of two bosons interacting via contact inter-particle potential trapped 
in the ID symmetric double-well potential 

Vsymmix) = O.Sx^ - 2.5V0.0784 + x2 + 3.16204 + H(A = 1, D = 0.75). (51) 

To obtain this potential we take the lowest-energy eigenvalue of Eq. ljH^ with Vi{x) = 0.5(x + 
2.5)2, V2ix) = 0.5(x - 2.5)2, C = 0.7, B = 0.0 and add a barrier function Vb{A = 1,D = 
0.75). The minimal value of this potential energy has been adjusted to zero by applying a 
constant shift. 

In FiglHlwe plot the ground state energy of this system as a function of the inter-particle 
interaction strength Aq, obtained within the framework of the self-consistent many-body 
MCHB(M) method, M=2,---,10. We also present the energies obtained within one- and 
two-orbital mean-fields, MF(1) and MF(2) respectively. All the energies are plotted with 
respect to the ground state energy E{0) of the non-interacting two-boson system. At this 
limit the ground state many-body wave-function is given by a single configuration where 
both bosons reside in the lowest-energy orbital of the respective bare Hamiltonian, and the 
total energy of this ground state E{0) is twice the respective orbital energy. 

In this figure it is very difficult to distinguish between the energy curves obtained within 
the two-orbital MCHB (2) theory and the multi-orbital MCHB(M) M=3,- ■ • ,10 ones. This 
observation tells us that for our example an adequate description can already be obtained 
within the two-orbital MCHB(2) method, the inclusion of more orbitals does not lead to 
significant "observable" improvements. We demonstrate this in the inset of FiglHl where 
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we present a part of the MCHB(M) energy curves on an enlarged scale. Comparing the 
energy gaps between successive MCHB(M) energy curves, we observe the convergence of the 
MCHB(M) method. 

To arrive at a deeper insight into the role of many-body effects we also study the trapped 
two-bosonic system within the framework of multi-orbital mean-field theory which is a lim- 
iting one-permanent case of the MCHB approach as we have shown in Sec ill CI We recall 
that the one-orbital mean-field MF(1) is the famous Gross-Pitaevskii mean-field. In Fig|Hl 
the MF(1) energy curve is depicted by a dashed line. The one-orbital mean-field solution 
describes a "condensate" where both bosons reside in the same orbital. The two-orbital 
mean-field MF(2) solution describes a situation where one boson resides in one orbital and 
the second boson occupies another orbital. For the two-boson system such a state can be 
considered as a two-fold "fragmented" state. Here we observe a "critical phenomenon". 
Up to some critical value of the inter-particle interaction strength the " condensed" solution 
is energetically more favorable than the "fragmented" one. From this critical Aq on, the 
ground state becomes two-fold "fragmented". The intersection of the MF(1) and MF(2) 
energy curves gives the critical interaction strength at Aq = 0.0203. In FiglHlwe mark this 
point by a cross. 

The MCHB theory gives the following many-body picture of this transition. The natural 
orbital analysis applied to the MCHB(M) solutions at each Aq reveals that the character of 
the many-body MCHB(M) ground state smoothly develops with inter-particle interaction 
strength from "condensed", where only one natural orbital is occupied, to the two-fold "frag- 
mented" where two natural orbital have dominant and nearly equal occupations. Having at 
hand the natural orbital occupation numbers as a function of Aq, we find the inflection point 
at Ao = 0.0142 and attribute it to the transition point. In FiglHlthis point is marked by 
a bold filled circle. The comparison between many-body and mean-field predictions shows 
that in contrast to the sharp transition obtained within the multi-orbital mean-fields, an 
inclusion of the many-body effect makes this transition smooth, i.e., it is a crossover. 

This investigation of the minimal many-body system has verified that qualitative pre- 
dictions on the transition from condensation to fragmentation in symmetric double wells 
can alr eady be obtained within the framework of the self-consistent multi-orbital mean-field 
theory j^H;!^- ^^^'^ demonstrated that the two-orbital MCHB(2) theory provides 
in this case an accurate quantitative description and the inclusion of more orbitals leads to 
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minor changes. Generally, MCHB(M) opens the door to treat any two-boson system. 

V. DISCUSSION AND CONCLUSIONS 

In this work we developed a complete variational many-body theory for systems of N 
trapped bosons interacting via a general two-body interaction potential. The many-body 
wave function of this system is expanded over orthogonal many-body basis functions (config- 
urations) . Each basis function is constructed as a symmetrized Hartree product (permanent) 
with N bosons distributed over M one-particle functions. These one-particle functions and 
the respective expansion coefficients are treated as the variational parameters in this the- 
ory. The optimal variational parameters are obtained self- consistently by solving a coupled 
system of non-eigenvalue intcgro-differential equations to get the one-particle functions and, 
by diagonalizing the secular Hamiltonian matrix problem, to find the expansion coefficients. 
To construct this matrix we derived general rules for matrix elements, which are of relevance 
also for other many-body theories. We call this self-consistent theory multi-configurational 
Hartree for bosons, or MCHB(M) where M specifies the number of one-particle functions 
involved. 

The properties of the MCHB(M) equations were discussed. These equations are formally 
also valid for any size of the many-body expansion, i.e., for any number of configurations 
used in the expansion. Therefore, the MCHB(M) theory allows to find also the best possible 
many-body solution within any restricted configurational space used. We have shown that in 
the limiting case where only one permanent forms the many-body expansion, the MCHB(M) 
theory boils down to the self-consistent mean fields. In the simplest case when all bosons 
reside in the same orbital one gets the Gross-Pitacvskii equation. The multi-orbital mean- 
field theory is obtained in the more general single-permanent case, when bosons are allowed 
to reside in several one-particle functions. 

We have shown that if the many-body basis set spans a complete subspace of the Hilbert 
space, namely, when all possible configurations appearing as permutations of N bosons over 
M orbitals are taken into account, then the MCHB(M) solution is invariant with respect to 
a unitary transformation (linear combination) of the MCHB(M) orbitals. This property has 
been used to demonstrate that eigenfunctions of the reduced one-particle density, i.e., the 
natural orbitals are the MCHB(M) solution as well. We proposed to analyze the ground 
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and excited states in the terms of the natural orbitals and natural occupation numbers to 
more easily identify the depletion and fragmentation of the condensates. 

In the second part of our work we implemented the MCHB(M) method with M—2 or- 
bitals. We applied it to study the ground and excited states of the bosonic systems with the 
popular contact inter-particle interaction trapped in one- and two-dimensional symmetric 
and asymmetric double- well traps. The considered configurational space was spanned by 
all possible permutations of N=1000 bosons over two orbitals. Two lowest-energy eigen- 
functions of the respective bare Hamiltonian were used to construct the often employed 
fixed-orbital many-body basis set. We compare the fixed-orbital many-body predictions 
with those obtained self-consistently via the MCHB(2) theory to investigate the impact and 
relevance of self-consistency. 

We performed several ground state studies of the bosonic system trapped in the sym- 
metric double-well trap. In the first study, we keep the shape of the symmetric double-well 
trap potential fixed and vary the strength Aq of the inter-particle interaction in order to 
study the ground state fragmentation. We have seen that self-consistent MCHB(2) theory 
predicts a gradual enhancement of the fragmented fraction with Aq up to some critical inter- 
particle interaction strength, where the fragmentation achieves its maximal value. Further 
increase of Aq causes gradual decreasing of the fragmentation, because the energy of the 
bosonic system in this regime becomes larger than the potential barrier. The many-body 
result obtained within fixed bare Hamiltonian functions predicts a gradual enhancement 
of the fragmentation with increasing Aq followed by an unphysical saturation of the frag- 
mented fraction to some constant value. In the second study, to investigate the transition 
point from condensation to fragmentation we keep the inter-particle interaction strength 
fixed and rump the barrier up. The critical value of the barrier height obtained with bare 
Hamiltonian functions are considerably underestimated in comparison with the more exact, 
self-consistent, MCHB(2) many-body predictions. A main conclusion derived from this in- 
vestigation is that the quantitative characterization of the ground state properties of the 
bosonic system trapped in symmetric double-wells can be obtained in the framework of 
self-consistent methods, the fixed-orbital many-body theory utilizing the same size of CI 
expansion can be unreliable even concerning qualitative predictions. 

We addressed the ground state properties of the bosonic system trapped in the asymmetric 
double-well potential. In this study we keep the shape of the asymmetric double-well trap 
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potential fixed and vary the strength of the inter-particle interaction. The MCHB(2) theory 
predicts two regimes for the ground state. In the first one the atomic cloud is localized in 
the deeper well; from some critical inter-particle interaction strength on the system enters 
the second regime where bosons occupy both wells. We show that such a picture can not 
be obtained within a fixed two-orbital many-body treatment. To overcome this difficulty, 
one must use at least three fixed orbitals to construct the permanents. However the active 
space, i.e., the number of many-body basis functions used in this case is substantial and 
often beyond reach. 

To exploit the freedom of unitary transformations we analyzed the distribution of the 
MCHB(2) expansion coefficients obtained for different linear combinations of the MCHB(2) 
orbitals for the ground state of the bosonic systems trapped in the symmetric and asym- 
metric double-wells. We verified that statistical means and variances can indeed be used to 
characterize adequately the distributions of the expansion coefficients. Moreover, we have 
seen that the distributions with minimal width, obtained by minimizing the variance, exhibit 
very smooth profiles, irrespective of the symmetry of the trap potential used. For the studied 
examples the profiles of the smooth minimal distributions can very well be approximated 
by continuous Gaussian functions. 

We also investigated the first excited state of the system trapped in the symmetric double- 
well and demonstrate that self-consistency can be very important for an adequate description 
of excited states. We used the natural orbitals analysis to classify the ground and excited 
states. It was shown that without self-consistency the lowest-energy excited state describes a 
depletion of the condensate. The inclusion of self-consistency significantly lowers the energy 
of the first excited state and on top of that drastically changes its character: Instead of 
describing a condensate with a slightly larger depleted fraction, it describes a fragmented 
condensate with a substantial degree of almost 30% fragmentation. 

As an illustrative example we investigate the ground state of a two dimensional bosonic 
system trapped in a symmetric double-well potential. We show that self-consistency is 
expected to be of high relevance for many-body studies on trapped bosonic systems also in 
higher dimensions. 

Finally, wc have shown that the two-orbital MCHB(2) theory can provide quite accurate 
quantitative description for the ground state of two-bosonic systems in symmetric double- 
well traps. The MCHB(M) has been implemented for two bosons and in an illustrative 
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example the inclusion of more orbitals leads only to minor changes. 



APPENDIX A: OFF-DIAGONAL ELEMENTS OF TWO-BODY DENSITY MA- 
TRIX 

In this appendix we evaluate the off-diagonal elements of the two-body density matrix 
Pijkl = {'^\blbjbkbi\'^). We use the same shorthand notations as defined in Sec lllB l1 where 
only occupation numbers of the involved orbitals are shown. For instance, the configuration 
(; Hi — 2; rij + 2; ) differs from n = {; nf, nj] ) by excitation of two bosons from (pi to (pj. In 
all cases different indices can not have the same value. 

Piijj = ^ClC(,n,.2;n,+2;)\ln,{n, - l){nj + l){nj + 2), 
n 

Piijk = ^QC(;„^_2;„^.+i;„,+i;)^nj(ni - 1 ) (tIj + l)(nfc + 1), 

n 

pijkl = ^QC(;„,^_l;„^_l;„, + l;„, + l;)Wninj-(nfc + + l), 



Pijkk 

= CHC(;n,-l:n,-l;n^+2;)yninj{nk + l)(nfc + 2), 



Piiij = ^QC(.„^_i;„^+i;)(ni - l)Jni{nj + 1), 



Pijjj = '^CtC(,ni-l;nj+l;)nj Jriiinj + 1), 



Pikkj — ''^C*iC(^.ni-l;nj+l;)nky ni{nj + 1). 
n 

All other matrix elements are zero or can be reduced to the above ones due to symmetries 
of the two-body operator: 

Pijkl = Pjikl = Pijlk = Pjilk- 

When the many-body function |\E') and one-particle orbitals are real functions, some addi- 
tional symmetries are implied: pu = P21, P2111 = P1112 P2221 = P1222 and P2211 = Pii22- 
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FIG. 1: (Color online) Right panel: The many-body energies per particle of the system of N=1000 
trapped bosons as a function of inter-particle interaction strength Aq obtained within many-body 
basis sets constructed on fixed (black) and self-consistent (red) orbitals. As the fixed orbitals we 
use the eigenfunctions of the bare Hamiltonian (BH) with potential plotted on left panel. The 
self-consistent orbitals have been obtained within the framework of the MCHB(2) method. To 
demonstrate better the impact of self-consistency on the ground state energy we plot as A the 
difference between both energy curves in the inset. Left panel: Geometry of the symmetric one- 
dimensional trap used, see Ea. H4U|) . It can be viewed as a combination of two harmonic traps 
separated by a barrier, see text for details. 
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FIG. 2: (Color online) Demonstration that a lack of self-consistency can have a drastic impact 
on predicted many-body properties of the ground state. Fragmentation is plotted as a function 
of the inter-particle interaction strength. Shown is the fragmentation of = 1000 bosons in the 
symmetric double-well potential of FigQ The occupation number p2 = N — pi oi the second 
natural orbital of the reduced one-particle density is plotted as a function of Aq. The solid (red) 
line with filled circles and the solid (black) line with open circles mark the self-consistent (MCHB) 
and fixed-orbital (BH) results, respectively. 
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FIG. 3: (Color online) Fragmentation as a function of barrier height. The occupation number 
P2 = N — pi the second natural orbital of the reduced one-particle density is plotted as a 
function of the barrier height for N=1000 bosons trapped in a symmetric double- well potential (see 
text). The solid (red) line with filled circles and the solid (black) line with open circles represent 
the self-consistent (MCHB) and fixed-orbital (BH) many-body results, respectively. To emphasize 
the "non-trivial" difference between self-consistent and fixed-orbital results, the difference between 
both curves is plotted as a dashed (blue) line. The level of ~ 25% of fragmentation is achieved at 
barrier height of Vsymm{^) ~ 6.5 units with fixed-orbital many-body method while self-consistency 
shifts this point to a barrier height Vsymm{^) ~ 10.5 units, see text for discussion. 
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FIG. 4: (Color online) Study of N=1000 bosons in an asymmetric double-well using the MCHB(2) 
method. The asymmetric double-well potential and three eigenfunctions of the respective bare 
Hamiltonians lowest in energy are schematically shown in the left panel. Depending on the inter- 
particle interaction strength Aq the ground state of the system can enter two different regimes. For 
A < Aq'~ the self-consistent MCHB(2) orbitals (solid lines) and respective normalized density (solid 
line with filled circles) are localized in the deeper well as depicted in the right lower panel. In the 
right upper panel it is shown that for Aq > Aq** the MCHB(2) natural orbitals are distributed over 
both wells. 
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FIG. 5: (Color online) The distribution of the ground state expansion coefficients C{ni,n2) obtained 
in MCHB(2) depends on the particular choice of the one-particle basis functions. The left panels 
show the results for N=1000 bosons trapped in a symmetric double-well (Aq = 0.01). The right 
panels show the analogous results for the bosons in an asymmetric double- well (Aq = 0.01). Lower 
panels refer to the case where natural orbitals are used to construct the many-body basis functions 
|ni,n2 >. The width of a distribution is characterized in terms of its variance o"^^ = {n\) — {ni)'^. By 
applying unitary transformations (rotations) on the natural orbitals, the width of the distributions 
of the expansion coefficients can be minimized. Upper panels show the obtained distributions of 
the expansion coefficients with the minimal widths. The minimal distributions in these systems are 



well approximated by continuous Gaussian functions 
solid lines, see text for details. 
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FIG. 6: (Color online) Demonstration that self-consistency can have a great impact on predicted 
many-body properties of excited states. Shown are results for N=1000 bosons trapped in a sym- 
metric double well potential (see text for details). Left panel: Energy levels of the ground and 
first excited states obtained self-consistently are labelled as "SC Cround State" and "SC Excited 
State", respectively. For comparison we present the energy level of the non-self-consistent first 
excited state labelled as "Excited state", obtained by using the optimized ground state orbitals. 
The occupation of the second natural orbital p2 {pi = N — P2) is indicated for each state. Right 
panel: the MCHB solutions (natural orbitals) of the self-consistent ground and excited states (solid 
lines). The respective normalized densities are plotted by solid lines with filled circles. 
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FIG. 7: (Color online) Illustration of the relevance of self-consistency for many-body treatments 
of bosonic systems in two-dimensions. Ground state density of N=1000 bosons at Aq = 0.01 in 
the symmetric two-dimensional double- well trap of Eq. (|5()|) . Due to the perfect symmetry of the 
trap potential, the total density consists of two equivalent parts each localized in one well. To 
highlight the difference, the part of the self-consistent MCHB(2) ground state density localized in 
the left well is plotted together with the part of the total density obtained with the fixed-orbital 
many-body method localized in the right well. 
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FIG. 8: (Color online) The total ground state energy of two bosons trapped in a symmetric double- 
well potential as a function of the inter-particle interaction strength Aq- All energies are plotted with 
respect to the ground state energy E{0) of the non-interacting system. The energy curves obtained 
within the framework of the two-orbital MCHB(2) and multi-orbital MCHB(M) M=3,- • • ,10 are 
very close to each other. To distinguish between these curves we enlarge the scale 100 times in 
the inset. To emphasize the role of many-body effects we also show the results obtained using 
the one-orbital MF(1) (Gross-Pitaevskii) and two-orbital MF(2) mean-fields. The cross marks 
the transition point between "condensation" and "fragmentation" defined as the intersection of 
the MF(1) and MF(2) energy curves. At the MCHB(M) level we observe a smooth development 
from "condensation" to "fragmentation" instead of the sharp transition; the filled circle marks the 
respective transition point. 
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